rm(list = ls())
library(foreign)
library(plyr)
library(sandwich)
library(data.table)
library(lmtest)

setwd("")

# Data comes from https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/EEP5MT
# See abortion_anon.tab
data <- read.dta('abortion_anon.dta')

#Abortion Results
data <- subset(data, data$respondent_t0==1)
# Compute factor analysis outcome
compute.factor.dv <- function(dv.names, respondent.booleans, print.loadings = TRUE){
  responders <- data[respondent.booleans,]
  # Factor analysis
  factor.obj <- princomp(responders[, dv.names], cor=TRUE)
  if(print.loadings) print(loadings(factor.obj))
  dv <- as.vector(factor.obj$scores[,1])
  # More positive values on the factor should indicate more tolerance; reverse otherwise.
  if(cor(dv, data$t0_abortion_factor[respondent.booleans], use="complete.obs") < 0) dv <- -1 * dv
  # Put in the order of the main data frame
  dv.in.order <- dv[match(data$anon_id, responders$anon_id)]
  # Rescale to mean 0 sd 1 in placebo group; treatment effects can then be interpreted
  # as the effect in standard deviations the treatment would have among an untreated
  # population.
  dv.in.order <- (dv.in.order - mean(dv.in.order[!data$treat_ind], na.rm=TRUE)) /
    sd(dv.in.order[!data$treat_ind], na.rm=TRUE)
  return(as.vector(dv.in.order))
}
# Compute the DVs in line with the above procedures.
all.dv.names.t1 <- c('t1_ballot_abortion', 't1_gssifthewomanshealthisse', 
                     't1_gssthewomanwantsitforany', 't1_gssifshebecamepregnantas', 
                     't1_gssifthereisastrongchanc', 't1_gssifthefamilyhasaverylo', 
                     't1_gssifsheisnotmarriedandd', 't1_gssifsheismarriedanddoes', 
                     't1_abort_itemitswomensfa', 't1_abort_itemwomenwhohav', 
                     't1_abort_itemwithmodernb', 't1_abort_itemtheresnothi', 
                     't1_abort_itemifawomanfee', 't1_abort_itemthegovernme', 
                     't1_abort_itemifafriendch')

# Omnibus DV of all primary outcomes.
data$all.dvs.t1 <- compute.factor.dv(all.dv.names.t1, data$respondent_t1==1 & !is.na(data$respondent_t1))
data <- subset(data, respondent_t1 == 1)
write.csv(data, "bks_data.csv")